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Abstract 

The effect of far-held boundary conditions on the evolution of a hnite- 
amplitude two-dimensional wave in the Blasius boundary layer is assessed. 
With the use of the parabolized stability equations (PSE) theory for the 
numerical computations, either asymptotic, Dirichlet, Neumann or mixed 
boundary conditions are imposed at various distances from the wall. The re- 
sults indicate that asymptotic and mixed boundary conditions yield the most 
accurate mean-how distortion and unsteady instability modes in comparison 
with the results obtained with either Dirichlet or Neumann conditions. 
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1. INTRODUCTION 


In a direct numerical simulation (DNS) of spatially growing disturbances 
in boundary-layer flows, the infinite domain in the streamwise direction x 
must be truncated to a finite length, while the semi-infinite extent in the 
plate-normal direction y may or may not be truncated. In past simulations, 
either a truncated domain in y [1-3] has been used or the semi-infinite domain 
has been mapped to a finite one with a change of variables [4-7]. 

The physical problem unequivocally prescribes the boundary conditions 
at y — > oo. Asymptotic boundary conditions can be implemented when the 
domain is truncated at a location y where the mean-flow has reached a con- 
stant value. At higher y locations the asymptotic boundary conditions can 
be approximated rather well by mixed boundary conditions (i.e. involving 
the function and its derivative and sometimes called Robin conditions). Sim- 
pler conditions such as homogeneous Neumann, or Dirichlet conditions can 
also be used; however, strictly speaking, these conditions are incorrect. The 
question is whether the loss of accuracy due to these more easily implemented 
but approximate boundary conditions is acceptable. 

Here, we look at the effect of the far-field boundary conditions on the 
evolution of a finite-amplitude two-dimensional wave in the Blasius bound- 
ary layer. We select either asymptotic, mixed, homogeneous Dirichlet, or 
homogeneous Neumann conditions and impose these conditions at various 
distances y max from the wall. For this study, we employ the parabolized 
stability equations (PSE) to take advantage of the low computational cost 
(e.g., typically 15 min on a workstation for a fully nonlinear two-dimensional 
calculation). 
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The frequency, the streamwise starting location xo, and the initial ampli- 
tude of the Tollmien-Schlichting (TS) wave used here were previously used 
by Bertolotti, Herbert, and Spalart [8] in a comparison of PSE and DNS 
results. The DNS code of Spalart mapped the infinite domain in y to a fi- 
nite domain via a mapping, which enforced the correct boundary condition 
at infinity and avoided the approximations caused by a truncated domain. 
Comparisons of results between PSE and DNS showed excellent agreement 
for all modes, including the mean-flow distortion [9]. 

In the comparison between DNS and PSE results by Joslin, Streett, and 
Chang [10], a discrepancy was found in the mean-flow distortion component. 
This discrepancy was attributed to the different far-field boundary condi- 
tions imposed in the two codes; the DNS code used homogeneous Dirichlet 
conditions, and the PSE code used homogeneous Dirichlet conditions for 
all unsteady modes and a homogeneous Neumann condition for the steady 
mean-flow distortion term. This discrepancy motivated the current study. 

2. GOVERNING EQUATIONS 

The reference length is 8(x o) = \J vx$ /U 0 Q , which is defined at the 
streamwise location xq. The corresponding Reynolds number at xo is Ro = 
UqoS(xo) /v = 400. The nondimensional frequency of the two-dimensional 
TS wave is F = 2 x 10 6 7 rfis/U^ = 86, which yields u = 0.0344. 

The PSE equations used in this work are described in references 8 and 
9. The incompressible disturbance equations are reduced to two variables 
(the u and v components of velocity) by taking the curl of the Navier- 
Stokes equations, eliminating the w velocity component with the conti- 
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nuity equation, and reducing the three governing equations to two equa- 
tions. The disturbance field is expanded in a six-term Fourier series in time, 
0,ut,2ut, ...,k cut, ...,5 cut, which leads to complex Fourier coefficients u k and 
D k for the velocity components. These equations take the form, 

[L+^N]q k + M^ = i? k k = 0, 1, 2, ..., 5 (1) 

ax ax 

where q k = {D k ,D k } is the vector of profile functions, a k is the complex 
wavenumber for mode k composed of a real part y k describing the growth 
rate and an imaginary part k a describing the wavenumber, the operators 
L,M,N depend on a k and frequency kuj, and contain derivatives only in y. 
The operator L contains the Orr-Sommerfeld and Squire operators, which 
are well known in the parallel-flow stability theory. The right-hand-side term 
i? k is the convolution term stemming from the nonlinear products. 

Introducing the finite difference form dq^/dx — > (qk - q ± ld )/d x and 
da^/dx — > (a k — a^ ld )/dx into equation (1) yields 

old 1 1 

[L + k k N + — M] q k = Rk + -Mqf. k = 0, 1, 2, ..., 5 (2) 

ax ax ax 

These coupled set of equations can be solved by marching in x. 

The physical boundary conditions at y — > oo impose vanishing D k and D k 
velocities for the unsteady modes (k > 0) and vanishing u o and constant fo 
velocities for the steady mode (k = 0). A finite fo allows for changes in the 
displacement thickness as the flow transitions to turbulence. This condition 
becomes Do = 0 in flows over bodies with a bounded streamwise extent. 

Outside the boundary layer, the operators L,M,N have constant co- 
efficients, and the solution decays exponentially. When the computational 
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domain is truncated in this region, boundary conditions can be imposed that 
yield the exact solution at all interior points. These boundary conditions 
have been presented by Keller [11] for the Orr-Sommerfeld operator, and are 
extended here to the PSE equation. The basic idea is to require that the so- 
lution has a null projection onto the subspace spanned by the exponentially 
growing eigensolutions of the operator on the left-hand side of equation ( 2 ). 
These eigensolutions are evaluated when equation (2) is re-written as a first 
order system. In the present formulation the highest derivative in y of u is 
2, and of v is 4, hence we introduce the vector Xk = {hk, Dk, , v^' } 
were the prime denotes differentiation w.r.t. y, and re-write ( 2 ) as 

A x k + B = r k (y) (3) 

dy 

where the matrices A and B depend on k and contain, in addition to the 
information in equation ( 2 ), the relations d(uk)/dy = d(vk)/dy = 0 ^, 
d(v^)/dy = 0^, and d(v^)/dy = v£' . We then compute, for each mode k, 
the eigensolutions {Ape;} that solve [A T — A;B T ]ei = 0, where T denotes 
transpose, and relabel these eigensolutions so that Ai,A 2 and A 3 have a 
negative real part. The requirement of zero projection onto the growing 
eigenmodes yields the following asymptotic boundary conditions, 


x k • B T ei = c k • e; i = 1,2,3 (4) 

where 

_ r h b 2 

Ct ~ { (c, +A,)’ (C2 + Ai)’''' } (5) 

is obtained by approximating each component of the forcing rk in the neigh- 
borhood of the boundary by a function of the form b e cy , and the dot product 
in equation (4) is defined as a • b = Y] a i h ■ 
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The spatial DNS code solves the disturbance form of the full Navier- 
Stokes equations with high-order finite- and compact-difference methods and 
spectral methods. Homogeneous Dirichlet conditions are imposed in the far- 
held and at the wall, inflow conditions consist of the Blasius and eigenfunc- 
tions provided by linear stability theory, and the buffer domain technique 
[12] is used for the streamwise outflow condition. Refer to references 13 for a 
discussion of accuracy issues with grid refinement and outflow buffer domain 
treatment. 


3. RESULTS 

Results were obtained for cases with the upper boundary placed at 
j/max = 15, 20, 30, 45, 60, 90, and 130. For each of these cases, com- 
putations were made that employed the following conditions for modes 
k = 0,1, 2, 3, 4, 5: 

Dirichlet conditions: 


Vk = 0 , v k = 0 , 


dv k 

dy 


0 


Neumann conditions: 


diik _ „ dv\ _ 

dy ~ ’ dy ~ ’ 


d 2 v k 
dy 2 


0 


(6) 


(7) 


Mixed conditions: 


duk 

dy 


+ dkllk 


0, 


dv k 
dy 


+ a k Vk 


0, 



+ a k 


dy k 
dy 


0 


(8) 


Asymptotic conditions: 


rri 

Xk • B e; = Ck • e; i = 1, 2, 3 


(9) 
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The conditions for the highest derivative of v with respect to y are 
derived from the continuity equation. The mixed boundary conditions and 
the asymptotic boundary conditions are altered for the mean-flow distortion 
term (i.e. k = 0), to the form discussed above, namely, 


u 0 = 0 , 


dv 0 

dy 


0, 


d 2 v 0 
dy 2 


0 


( 10 ) 


The initial condition was composed of the single Fourier mode k = 1, 
with an amplitude of 0.25 percent rms. based on the maximum of the u 
component of velocity. 

High resolution in the plate-normal direction y was obtained with five 
subdomains. In each subdomain, the u k and 0k velocity components were 
expanded with 18 Chebyshev polynomials. The step size in x was set to 
Ax = 10. 

Figure 1 shows the evolution of the disturbance amplitude based on the 
u component of velocity for the Fourier modes F = 1, F = 2, and the steady 
component F = 0 with a Reynolds number of R = U 00 6(x)lv = \JxR$ for 
results calculated by both the PSE and DNS codes. Both codes enforced 
the Dirichlet boundary conditions (eqs. 1) at y max = 130. The results agree 
well, which reasonably indicates the equivalence of the two procedures for 
the flat-plate problem. 

Next, computations were conducted with PSE theory to compare the 
maximum amplitudes of F = 0 and F = 1 modes as function of far-field 
boundary locations. At the downstream location that corresponds to R = 
940, Figure 2 displays the dependence of the maximum amplitudes, based on 
the u component of velocity, with the truncated far-field boundary location 
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J/max for the F = 1 and F = 0 modes. The solid line represents results that 
were obtained with asymptotic boundary conditions (eqs. (9)), the dashed 
line represents the mixed conditions (eqs. (8)), the square symbols represent 
Dirichlet conditions (eqs. (6)), the triangular symbols represent Neumann 
conditions (eqs. (7)), and the arrows denote the results obtained by applying 
the physical boundary conditions at infinity (using an algebraic mapping). 
Note that the boundary-layer edge (99% definition) grows from y = 5 at 
R = 400 to y = 12 at R = 940; therefore, the far-field boundary must be 
beyond y = 12. 

The results obtained with the asymptotic boundary conditions are in- 
dependent of y max once the mean-flow has reached a constant value (e.g. 
with less than a 0.01% variation). To increase accuracy, the operators in 
equation (3) can be evaluated with the mean-flow value at infinity, rather 
than at y max , however, for y max < 20 the mean-flow is still varying when the 
boundary is reached, and, consequently, the asymptotic boundary conditions 
become only approximate. The dip in the solid curve in figure 2 displays this 
fact. The evaluation of Ck in equation (4) can also affect accuracy; when the 
exponential fit (which is exact in the linear case) of rk is replaced by a two 
term Taylor series approximation of rk at y max , the calculated growth rates 
fall 15% short of the exact value, even for y max > 20. 

The mixed boundary conditions impose the exponential decay exp (—ay) 
to the solution. The complex wavenumber a is an eigenvalue of A T — A;B T , 
and when the other two decaying eigensolutions have a decay rate much 
higher than a, the mixed boundary conditions become equivalent to the 
asymptotic boundary conditions, provided y max is sufficiently large. A dif- 


7 



ference between the dashed and solid curves in figure 2 is barely visible at 
j/max < 35, although such a close agreement should not always be expected. 
However, for y max > 35, the mixed and asymptotic boundary conditions lead 
to the same solutions. Because the mixed conditions are homogeneous and 
have a simpler form, they are easier to implement. 

The Neumann conditions yield accurate results for y max > 45. These 
conditions allow for a change in the boundary-layer displacement thickness 
(i.e. , nonzero v velocity for F = 0 at y max ). Similarly, Dirichlet bound- 
ary conditions lead to accurate results for the traveling mode F = 1 when 
ymax > 45. However, the steady component F = 0 is adversely affected 
by the Dirichlet boundary conditions even for large values of y max , as indi- 
cated by the square symbols in Figure 2a. Furthermore, the v component 
of velocity vanishes, which prevents changes from the laminar value of the 
boundary-layer displacement thickness. Figure 3 displays the u$ and Vo ve- 
locity components for the mean-flow distortion mode (F = 0) at R = 940. 
Outside the boundary layer, the Vo velocity decreases linearly to match the 
zero boundary value at y max = 130. Because the boundary-layer displace- 
ment thickness tends to increase beyond the laminar value, mass conservation 
forces a nonzero u o component of velocity outside the boundary layer, which, 
in turn, creates an artificial boundary layer at y max . (Note that the small 
errors between the DNS and PSE v 0 profiles are of the order 10 — 6 , which 
can be attributed to numerical errors in the DNS approach.) 

All PSE calculations up to this point have been done using of 90 Cheby- 
shev polynomials in y per variable, per mode. This high resolution was 
chosen in order to remove resolution issues from the analysis. To assest the 



effect of the far-held resolution, we have repeated the computation using two 
domains, the inner one going from the wall to y = 5 and the outer one from 
y = 5 to y = 30. A linear mapping from physical to [—1, 1] was used in both 
domains. The resolution in the lower domain was fixed at 18 polynomials per 
variable, and in the outer domain either 10 or 30 polynomials were employed. 
Table 1 below displays the maximum u amplitude for the F = 0 and F = 1 
modes at R = 940 obtained with different boundary conditions. The exact 
values are 0.595 % for F = 0 and 2.843 % for F = 1. (The column labeled 
’’High Res” displays the values shown in figure 2.) 


Table 1. Modal maximums at R = 940 




O 

II 



F=1 


BC Type 

10 

30 

High Res 

10 

30 

High Res 

Asymptotic 

0.605 

0.595 

0.596 

2.966 

2.840 

2.844 

Mixed 

0.597 

0.597 

0.598 

2.869 

2.853 

2.858 

Neumann 

0.666 

0.684 

0.684 

3.709 

3.807 

3.812 

Dirichlet 

0.152 

0.286 

0.298 

2.010 

1.855 

1.895 


The asymptotic, mixed, and Neumann conditions display only a small 
variation with change in resolution. The Dirichlet condition is more sensitive, 
due to the need to resolve the artificial boundary-later at the upper domain, 
shown in figure 3. The largest difference between results, thus, comes from 
the boundary condition implemented, rather than the resolution. 

4. CONCLUSIONS 

The use of a finite domain in y plus Dirichlet and Neumann boundary 
conditions eliminates some coding difficulties in direct Navier- Stokes simula- 
tion codes, but introduces errors. As in the case considered here, the errors 
are small when the truncation location y max is located well into the region 
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of exponential decay of the disturbances. An exception is the steady com- 
ponent F = 0, which does not decay in the free stream and for which the 
error introduced by the use of Dirichlet conditions does not vanish as y max 
is increased. A similar error also is expected for three-dimensional steady 
disturbances because they decay slowly (i.e., as in exp (— /3 2 y)) in the free 
stream. The errors introduced in the calculation of traveling modes by ei- 
ther Dirichlet or Neumann conditions, on the other hand, are negligible if 
a truncation location y max is chosen sufficiently far from the plate. In con- 
trast, asymptotic boundary conditions and mixed boundary conditions yield 
accurate results when imposed beyond the 99.99% definition of the boundary 
layer edge. The asymptotic conditions are exact, but require a significantly 
greater amount of coding to implement. 
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FIG. 1. Amplitudes of F = 0, F = 1, and F = 2 modes from PSE (solid 
line) and DNS (symbols) with Dirichlet boundary conditions. 

FIG. 2. Maximum amplitudes of F = 0 and F = 1 modes as function of y m 
for the case of asymptotic (solid line), mixed (dashed line), Neumann (trian- 
gles), and Dirichlet boundary conditions (squares), and physical boundary 
conditions at y max — > oo (arrow) at R = 940. 

FIG. 3. Velocity components u and v for the F = 0 mode of PSE theory 
with asymptotic (line) and Dirichlet (dashed) boundary conditions and DNS 
results (symbols) at R = 940. 
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FIG. 2. Maximum amplitudes of F = 0 and F = 1 modes as function of y m 
for the case of asymptotic (solid line), mixed (clashed line), Neumann (trian- 
gles), and Dirichlet boundary conditions (squares), and physical boundary 
conditions at y max — > oo (arrow) at R = 940. 
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F=0 


F=0 



FIG. 3. Velocity components u and v for the F = 0 mode of PSE theory 
with asymptotic (line) and Dirichlet (dashed) boundary conditions and DNS 
results (symbols) at R = 940. 
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